function [oModel] = model_T_update(iModel,T,dt)
%MODEL_UPDATE update state from t to t+dt.
%
% [State] = model_update(w,State0,model_settings) 
%
% INPUTS:
%   iModel: a structure of the model
%   w: deposition rate, layer/s
%   T: temperature, K
%   dt: time step.
%
% OUTPUTS:
%   State = a structure with the following fields:
%    h
%    rho
%    meanAlpha2
%    meanBeta2
%    meanR2
%    varAlpha2
%    varBeta2
%    varR2
%    meanM2
%
% SEE ALSO: model_id.m, model_ini.m

%%
% assert(isstruct(State0),'State0 must be a structure');
% assert(all(isfield(State0,{'h','rho','meanAlpha2','meanBeta2','varAlpha2','varBeta2','meanR2','varR2','meanM2'})));
% assert(isstruct(model_settings),'model_settings must be a structure');
% assert(all(isfield(model_settings,{'FitTo','FittingTimeRange','W','nu','sigma2','K','Tau','dt','dx','M2ModeWeighting','dt','dx','L0','L'})));

%%
oModel = iModel;

%% Linear interpolation to get model coefficient
% Tab_DepRate   = oModel.W;
% Tab_R_nu      = oModel.nu;
% Tab_R_sigma2  = oModel.sigma2;
% Tab_R_Rh      = oModel.Rh;
% Tab_D_K       = oModel.K;
% Tab_D_Tau     = oModel.Tau;

Tab_T         = oModel.T;
Tab_R_nu      = oModel.nu;
Tab_R_sigma2  = oModel.sigma2;
Tab_R_Rh      = oModel.Rh;
Tab_D_K       = oModel.K;
Tab_D_Tau     = oModel.Tau;

if (Tab_T(1) < T) && (T <= Tab_T(end))
    nu     = interp1(Tab_T,Tab_R_nu,T,'linear','extrap');
    sigma2 = interp1(Tab_T,Tab_R_sigma2,T,'linear','extrap');
    Rh     = interp1(Tab_T,Tab_R_Rh,T,'linear','extrap');
    K      = interp1(Tab_T,Tab_D_K,T,'linear','extrap');
    tau    = interp1(Tab_T,Tab_D_Tau,T,'linear','extrap');
elseif T <= Tab_T(1)
    nu     = Tab_R_nu(1);
    sigma2 = Tab_R_sigma2(1);
    Rh     = Tab_R_Rh(1);
    K      = Tab_D_K(1);
    tau    = Tab_D_Tau(1);
else
    nu     = Tab_R_nu(end);
    sigma2 = Tab_R_sigma2(end);
    Rh     = Tab_R_Rh(end);
    K      = Tab_D_K(end);
    tau    = Tab_D_Tau(end);
end

%%
% The following models of h, rho and roughness square (R2) are documented
% in detail in Hu2009_CDC

% Thickness ---------------------------------------------------------
oModel.h = iModel.h+Rh*dt;
mode = oModel.mode;
% SOR  --------------------------------------------------------------
oModel.rho = (iModel.rho*iModel.h+Rh*(K*dt+(K-iModel.rho)*tau*(exp(-dt/tau)-1)))/(iModel.h+Rh*dt);

% varAlpha2 and varBeta2 --------------------------------------------
% NOTE: Calculate varAlpha2 and varBeta2 first because meanAlpha2 and
% meanBeta2 will be updated.
oModel.varAlpha2 = zeros(mode,1);
oModel.varBeta2  = zeros(mode,1);
for j = 1:mode
    temp1 = exp(-2*nu*j^2*dt);
    temp2 = exp(-4*nu*j^2*dt);
    temp3 = sigma2*(temp1-1)/(-2*nu*j^2);
    oModel.varAlpha2(j) = temp2*iModel.varAlpha2(j)+4*temp1*temp3*iModel.meanAlpha2(j)+2*temp3^2;
    oModel.varBeta2(j)  = temp2*iModel.varBeta2(j)+4*temp1*temp3*iModel.meanBeta2(j)+2*temp3^2;
end
oModel.varR2 = sum(oModel.varAlpha2+oModel.varBeta2)/(4*pi^2);

% meanAlpha2 and meanBeta2 ------------------------------------------
oModel.meanAlpha2 = zeros(mode,1);
oModel.meanBeta2  = zeros(mode,1);
for j = 1:mode
    temp = sigma2/(2*nu*j^2);
    temp2 = exp(-2*nu*dt*j^2);
    oModel.meanAlpha2(j) = temp+(iModel.meanAlpha2(j)-temp)*temp2;
    oModel.meanBeta2(j)  = temp+(iModel.meanBeta2(j)-temp)*temp2;
end
oModel.meanR2 = model_CalMeanR2(oModel);
oModel.meanM2 = model_CalMeanM2(oModel);